COVID19 GER

Prepare data

Prevalence

df_ger_prev <- read_csv('Germany_timeseries_prep.csv')

Parsed with column specification:
cols(
  date = col_character(),
  anzahlfall = col_double(),
  kreis = col_double(),
  ewz = col_double(),
  shape__area = col_double(),
  cumcase = col_double(),
  kreis_name = col_character(),
  extra = col_double(),
  agree = col_double(),
  sci = col_double(),
  neuro = col_double(),
  open = col_double(),
  rate_day = col_double(),
  popdens = col_double(),
  runday = col_double()
)

df_ger_prev <- df_ger_prev %>% mutate(date = as.Date(date, "%d%b%Y"),
                                  kreis = as.character(kreis)) %>% 
  dplyr::select(kreis, date, rate_day)

df_ger_prev

Scoial distancing

df_ger_socdist <- read_csv('Germany_socdist_fb_kreis.csv')

Parsed with column specification:
cols(
  all_day_ratio_single_tile_users = col_double(),
  date = col_date(format = ""),
  kreis = col_double(),
  kreis_name = col_character(),
  extra = col_double(),
  agree = col_double(),
  sci = col_double(),
  neuro = col_double(),
  open = col_double(),
  frequ = col_double(),
  runday = col_double()
)

df_ger_socdist <- df_ger_socdist %>% 
  rename(socdist_single_tile = all_day_ratio_single_tile_users) %>%
  select(kreis, date, socdist_single_tile) %>% 
  mutate(kreis = as.character(kreis))

df_ger_socdist

NA

Personality

df_ger_pers <- read_csv('Germany_timeseries_prep.csv')

Parsed with column specification:
cols(
  date = col_character(),
  anzahlfall = col_double(),
  kreis = col_double(),
  ewz = col_double(),
  shape__area = col_double(),
  cumcase = col_double(),
  kreis_name = col_character(),
  extra = col_double(),
  agree = col_double(),
  sci = col_double(),
  neuro = col_double(),
  open = col_double(),
  rate_day = col_double(),
  popdens = col_double(),
  runday = col_double()
)

df_ger_pers <- df_ger_pers %>% 
  select(kreis, open, sci, extra, agree, neuro) %>%
  dplyr::rename(pers_o = open, 
         pers_c = sci,
         pers_e = extra,
         pers_a = agree,
         pers_n = neuro) %>%
  distinct() %>%
  mutate(kreis = as.character(kreis))

df_ger_pers

NA

Controls

df_ger_ctrl <- read.csv2('Germany_controls.csv', sep = ';', dec=',')

df_ger_ctrl <- df_ger_ctrl %>% select(-kreis_nme) %>%
    mutate(kreis = as.character(kreis),
           popdens = popdens %>% 
             as.character() %>%
             str_replace('\\.', '')%>%
             as.numeric())

df_ger_ctrl

NA

Merge prevalence data

# create sequence of dates
date_sequence <- seq.Date(min(df_ger_prev$date),
                          max(df_ger_prev$date), 1)
                     
# create data frame with time sequence
df_dates = tibble(date_sequence, 1:length(date_sequence)) 
names(df_dates) <- c('date', 'time')


# merge prevalence data
df_ger_prev <- df_ger_prev %>% 
  inner_join(df_ger_pers, by = 'kreis') %>%
  inner_join(df_ger_ctrl, by = 'kreis') %>%
  merge(df_dates, by='date') %>% 
  arrange(kreis)

df_ger_prev

Merge socdist data

# create sequence of dates
date_sequence <- seq.Date(min(df_ger_socdist$date),
                          max(df_ger_socdist$date), 1)
                     
# create data frame with time sequence
df_dates = tibble(date_sequence, 1:length(date_sequence)) 
names(df_dates) <- c('date', 'time')

# merge socdist data
df_ger_socdist <- df_ger_socdist %>% 
  inner_join(df_ger_pers, by = 'kreis') %>%
  inner_join(df_ger_ctrl, by = 'kreis') %>% 
  merge(df_dates, by='date') %>% 
  arrange(kreis)

df_ger_socdist

NA

Control for weekend effect

easter <- seq.Date(as.Date('2020-04-10'), as.Date('2020-04-13'), 1)

df_ger_loess <- df_ger_socdist %>% 
  mutate(weekday = format(date, '%u')) %>% 
  filter(!(weekday %in% c('6','7') | date %in% easter)) %>% 
  split(.$kreis) %>%
  map(~ loess(socdist_single_tile ~ time, data = .)) %>%
  map(predict, 1:max(df_ger_socdist$time)) %>% 
  bind_rows() %>% 
  gather(key = 'kreis', value = 'loess') %>% 
  group_by(kreis) %>% 
  mutate(time = row_number())

df_ger_socdist <- df_ger_socdist %>% merge(df_ger_loess, by=c('kreis', 'time')) %>% 
  mutate(weekday = format(date, '%u')) %>% 
  mutate(socdist_single_tile_clean = ifelse(weekday %in% c('6','7') | date %in% easter, 
                                            loess, socdist_single_tile)) %>%
  arrange(kreis, time) %>% 
  select(-weekday)

df_ger_socdist %>% drop_na()

NA

Explore data

Plot prevalence over time

df_ger_prev %>% ggplot(aes(x=time, y=rate_day)) + 
  geom_point(aes(col=kreis, size=popdens)) + 
  geom_smooth(method="loess", se=T) + 
  theme(legend.position="none") +
  ggtitle("Overall prevalence over time")

pers <- c('pers_o', 'pers_c', 'pers_e', 'pers_a', 'pers_n')

for (i in pers){

gg <- df_ger_prev %>% mutate(prev_tail = cut(.[[i]], 
      breaks = c(-Inf, quantile(.[[i]], 0.2), quantile(.[[i]], 0.8), Inf),
      labels = c('lower tail', 'center', 'upper tail'))) %>% 
  filter(prev_tail != 'center') %>%
  ggplot(aes(x=time, y=rate_day)) + 
  geom_point(aes(col=kreis, size=popdens)) + 
  geom_smooth(method="loess", se=T) + 
  facet_wrap(~prev_tail) + 
  theme(legend.position="none") +
  ggtitle(i)

print(gg)
}

Plot social distancing (single tile) over time

df_ger_socdist %>% ggplot(aes(x=time, y=socdist_single_tile_clean)) + 
  geom_point(aes(col=kreis, size=popdens)) + 
  geom_smooth(method="loess", se=T) + 
  theme(legend.position="none") +
  ggtitle("Overall social distancing (single tile) over time")

pers <- c('pers_o', 'pers_c', 'pers_e', 'pers_a', 'pers_n')

for (i in pers){

gg <- df_ger_socdist %>% mutate(prev_tail = cut(.[[i]], 
      breaks = c(-Inf, quantile(.[[i]], 0.2), quantile(.[[i]], 0.8), Inf),
      labels = c('lower tail', 'center', 'upper tail'))) %>% 
  filter(prev_tail != 'center') %>%
  ggplot(aes(x=time, y=socdist_single_tile_clean)) + 
  geom_point(aes(col=kreis, size=popdens)) + 
  geom_smooth(method="loess", se=T) + 
  facet_wrap(~prev_tail) + 
  theme(legend.position="none") +
  ggtitle(i)

print(gg)
}

df_ger_socdist <- df_ger_socdist %>% mutate(socdist_single_tile = socdist_single_tile_clean) %>% 
  select(-loess, -socdist_single_tile_clean)

Correlations

df_ger_prev %>% group_by(kreis) %>% 
  summarize_if(is.numeric, mean) %>% 
  select(-kreis, -time) %>% 
  cor(use='pairwise.complete') %>% 
  round(3) %>% as.data.frame()

df_ger_socdist %>% group_by(kreis) %>% 
  summarize_if(is.numeric, mean) %>% 
  select(-kreis, -time) %>% 
  cor(use='pairwise.complete') %>% 
  round(3) %>% as.data.frame()

NA
NA

Rescale data

lvl2_scaled <- df_ger_prev %>% 
  dplyr::select(-time, -rate_day, -date) %>% 
  distinct() %>% 
  mutate_at(vars(-kreis), scale)
  
lvl1_scaled <- df_ger_prev %>% select(kreis, time, rate_day)

df_ger_prev_scaled <- plyr::join(lvl1_scaled, lvl2_scaled, by = 'kreis')

df_ger_prev_scaled

NA

lvl2_scaled <- df_ger_socdist %>% 
  dplyr::select(-time, -socdist_single_tile, -date) %>% 
  distinct() %>% 
  mutate_at(vars(-kreis), scale)
  
lvl1_scaled <- df_ger_socdist %>% select(kreis, time, socdist_single_tile)

df_ger_socdist_scaled <- plyr::join(lvl1_scaled, lvl2_scaled, by = 'kreis')

df_ger_socdist_scaled

NA

Predict Prevalence

Extract first day of covid outbreak

# get onset day
df_ger_onset_prev <- df_ger_prev_scaled %>% 
  group_by(kreis) %>% 
  mutate(rate_cs = cumsum(rate_day)) %>% 
  filter(rate_cs > 0) %>%
  summarize(onset_prev = min(time))
  
# merge with county data
df_ger_onset_prev <- df_ger_prev_scaled %>% 
  select(-time, -rate_day) %>%
  distinct() %>% 
  left_join(df_ger_onset_prev, by = 'kreis')

# handle censored data
df_ger_onset_prev <- df_ger_onset_prev %>% 
  mutate(event = ifelse(is.na(onset_prev), 0, 1)) %>% 
  mutate(onset_prev = replace_na(onset_prev, as.numeric(diff(range(df_ger_prev$date)))+1))

df_ger_onset_prev

NA

Extract slopes

# cut time series before onset
df_ger_prev_scaled <- df_ger_prev_scaled %>% 
  group_by(kreis) %>% 
  mutate(rate_cs = cumsum(rate_day)) %>% 
  filter(rate_cs > 0) %>%
  mutate(time = time-min(time)+1) %>%
  ungroup() %>%
  filter(time <= 30) %>%
  select(-rate_cs)

# drop counties with little data
df_ger_prev_scaled <- df_ger_prev_scaled %>%
  group_by(kreis) %>%
  filter(n() == 30) %>%
  ungroup()

# extract slope prevalence
df_ger_slope_prev <- df_ger_prev_scaled %>% 
  split(.$kreis) %>% 
  map(~ lm(rate_day ~ time, data = .)) %>%
  map(coef) %>% 
  map_dbl('time') %>% 
  as.data.frame() %>% 
  rownames_to_column('kreis') %>% 
  rename(slope_prev = '.')

# merge with county data
df_ger_slope_prev <- df_ger_prev_scaled %>% 
  select(-time, -rate_day) %>%
  distinct() %>% 
  inner_join(df_ger_slope_prev, by = 'kreis') %>%
  drop_na()

Explore distributions

df_ger_onset_prev %>% ggplot(aes(onset_prev)) + geom_histogram()

df_ger_slope_prev %>% ggplot(aes(slope_prev)) + geom_histogram()

Predict COVID onset with time-to-event regression

# predict onset from personality
cox_onset_prev <- coxph(Surv(onset_prev, event) ~ 
                          pers_o + pers_c + pers_e + pers_a + pers_n, 
                        data = df_ger_onset_prev)
cox_onset_prev %>% summary()

Call:
coxph(formula = Surv(onset_prev, event) ~ pers_o + pers_c + pers_e + 
    pers_a + pers_n, data = df_ger_onset_prev)

  n= 400, number of events= 400 

            coef exp(coef)  se(coef)      z Pr(>|z|)    
pers_o  0.131653  1.140712  0.051953  2.534   0.0113 *  
pers_c -0.104500  0.900775  0.053314 -1.960   0.0500 *  
pers_e  0.056186  1.057794  0.052903  1.062   0.2882    
pers_a -0.008624  0.991413  0.049615 -0.174   0.8620    
pers_n -0.276799  0.758207  0.057728 -4.795 1.63e-06 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

       exp(coef) exp(-coef) lower .95 upper .95
pers_o    1.1407     0.8766    1.0303     1.263
pers_c    0.9008     1.1102    0.8114     1.000
pers_e    1.0578     0.9454    0.9536     1.173
pers_a    0.9914     1.0087    0.8995     1.093
pers_n    0.7582     1.3189    0.6771     0.849

Concordance= 0.592  (se = 0.018 )
Likelihood ratio test= 39.92  on 5 df,   p=2e-07
Wald test            = 37.99  on 5 df,   p=4e-07
Score (logrank) test = 37.99  on 5 df,   p=4e-07

# predict onset from personality with controls
cox_onset_prev_ctrl <- coxph(Surv(onset_prev, event) ~ 
                               pers_o + pers_c + pers_e + pers_a + pers_n +
                               women + academics + cdu + afd + hospital_beds +
                               tourism_beds + gdp + manufact + airport + age +
                               popdens,
                             data = df_ger_onset_prev)
cox_onset_prev_ctrl %>% summary()

Call:
coxph(formula = Surv(onset_prev, event) ~ pers_o + pers_c + pers_e + 
    pers_a + pers_n + women + academics + cdu + afd + hospital_beds + 
    tourism_beds + gdp + manufact + airport + age + popdens, 
    data = df_ger_onset_prev)

  n= 392, number of events= 392 
   (8 observations deleted due to missingness)

                   coef exp(coef)  se(coef)      z Pr(>|z|)    
pers_o         0.060631  1.062506  0.059458  1.020 0.307861    
pers_c        -0.082813  0.920523  0.060718 -1.364 0.172598    
pers_e         0.009496  1.009541  0.059390  0.160 0.872972    
pers_a        -0.019747  0.980447  0.057978 -0.341 0.733413    
pers_n        -0.249887  0.778889  0.064236 -3.890 0.000100 ***
women          0.083662  1.087261  0.068554  1.220 0.222322    
academics      0.335157  1.398160  0.100092  3.348 0.000813 ***
cdu            0.179656  1.196805  0.079129  2.270 0.023182 *  
afd           -0.017508  0.982644  0.079368 -0.221 0.825408    
hospital_beds -0.233707  0.791594  0.068731 -3.400 0.000673 ***
tourism_beds   0.055041  1.056584  0.055716  0.988 0.323208    
gdp           -0.151600  0.859332  0.115313 -1.315 0.188618    
manufact       0.112626  1.119214  0.095688  1.177 0.239190    
airport       -0.179572  0.835627  0.062056 -2.894 0.003807 ** 
age           -0.102691  0.902406  0.090956 -1.129 0.258890    
popdens        0.125177  1.133349  0.087247  1.435 0.151360    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

              exp(coef) exp(-coef) lower .95 upper .95
pers_o           1.0625     0.9412    0.9456    1.1938
pers_c           0.9205     1.0863    0.8172    1.0369
pers_e           1.0095     0.9905    0.8986    1.1342
pers_a           0.9804     1.0199    0.8751    1.0984
pers_n           0.7789     1.2839    0.6867    0.8834
women            1.0873     0.9197    0.9506    1.2436
academics        1.3982     0.7152    1.1491    1.7012
cdu              1.1968     0.8356    1.0249    1.3976
afd              0.9826     1.0177    0.8411    1.1480
hospital_beds    0.7916     1.2633    0.6918    0.9057
tourism_beds     1.0566     0.9464    0.9473    1.1785
gdp              0.8593     1.1637    0.6855    1.0772
manufact         1.1192     0.8935    0.9278    1.3501
airport          0.8356     1.1967    0.7399    0.9437
age              0.9024     1.1081    0.7551    1.0785
popdens          1.1333     0.8823    0.9552    1.3447

Concordance= 0.659  (se = 0.016 )
Likelihood ratio test= 106.5  on 16 df,   p=2e-15
Wald test            = 104.3  on 16 df,   p=5e-15
Score (logrank) test = 108.4  on 16 df,   p=9e-16

Predict prevalence slopes with linear models

# predict slopes from personality
lm_slope_prev <- lm(slope_prev ~ pers_o + pers_c + pers_e + pers_a + pers_n, 
                         data = df_ger_slope_prev)
lm_slope_prev %>% summary()

Call:
lm(formula = slope_prev ~ pers_o + pers_c + pers_e + pers_a + 
    pers_n, data = df_ger_slope_prev)

Residuals:
      Min        1Q    Median        3Q       Max 
-0.032829 -0.015670 -0.005169  0.008046  0.235669 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  0.0276072  0.0013833  19.957   <2e-16 ***
pers_o      -0.0008705  0.0014835  -0.587    0.558    
pers_c      -0.0016391  0.0015871  -1.033    0.302    
pers_e       0.0021386  0.0015724   1.360    0.175    
pers_a       0.0024992  0.0015982   1.564    0.119    
pers_n       0.0016993  0.0016343   1.040    0.299    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.02738 on 386 degrees of freedom
Multiple R-squared:  0.01215,   Adjusted R-squared:  -0.0006472 
F-statistic: 0.9494 on 5 and 386 DF,  p-value: 0.4489

# predict slopes from personality with controls
lm_slope_prev_ctrl <- lm(slope_prev ~  
                               pers_o + pers_c + pers_e + pers_a + pers_n +
                               women + academics + cdu + afd + hospital_beds +
                               tourism_beds + gdp + manufact + airport + age +
                               popdens,
                         data = df_ger_slope_prev)
lm_slope_prev_ctrl %>% summary()

Call:
lm(formula = slope_prev ~ pers_o + pers_c + pers_e + pers_a + 
    pers_n + women + academics + cdu + afd + hospital_beds + 
    tourism_beds + gdp + manufact + airport + age + popdens, 
    data = df_ger_slope_prev)

Residuals:
      Min        1Q    Median        3Q       Max 
-0.047064 -0.014053 -0.002229  0.008977  0.217515 

Coefficients:
                Estimate Std. Error t value Pr(>|t|)    
(Intercept)    0.0274783  0.0013022  21.101  < 2e-16 ***
pers_o        -0.0003747  0.0016078  -0.233   0.8159    
pers_c        -0.0001953  0.0015584  -0.125   0.9004    
pers_e         0.0013162  0.0015344   0.858   0.3915    
pers_a         0.0012651  0.0015783   0.802   0.4233    
pers_n         0.0008380  0.0015918   0.526   0.5989    
women          0.0046402  0.0019204   2.416   0.0162 *  
academics      0.0005925  0.0025926   0.229   0.8194    
cdu            0.0096133  0.0020618   4.663 4.35e-06 ***
afd            0.0026729  0.0020254   1.320   0.1877    
hospital_beds  0.0026756  0.0017517   1.527   0.1275    
tourism_beds  -0.0023895  0.0014941  -1.599   0.1106    
gdp            0.0005674  0.0030072   0.189   0.8504    
manufact       0.0028811  0.0024162   1.192   0.2339    
airport        0.0033653  0.0016388   2.054   0.0407 *  
age           -0.0022077  0.0025050  -0.881   0.3787    
popdens       -0.0003063  0.0022367  -0.137   0.8911    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.02575 on 375 degrees of freedom
Multiple R-squared:  0.1509,    Adjusted R-squared:  0.1146 
F-statistic: 4.164 on 16 and 375 DF,  p-value: 2.01e-07

CRF predicting slopes

ctrls <- cforest_unbiased(ntree=500, mtry=5)

crf_slope_prev <- cforest(slope_prev ~  
                               pers_o + pers_c + pers_e + pers_a + pers_n +
                               women + academics + cdu + afd + hospital_beds +
                               tourism_beds + gdp + manufact + airport + age +
                               popdens,
                           data = df_ger_slope_prev, 
                         controls = ctrls)

crf_slope_prev_varimp <- varimp(crf_slope_prev, nperm = 1)
crf_slope_prev_varimp_cond <- varimp(crf_slope_prev, conditional = T, nperm = 1)

crf_slope_prev_varimp %>% as.data.frame() %>% 
  rownames_to_column('variable') %>%
  ggplot(aes(x=variable, y=.)) +
  geom_bar(stat = 'identity') + 
  theme(axis.text.x = element_text(angle = 90))

crf_slope_prev_varimp_cond %>% as.data.frame() %>% 
  rownames_to_column('variable') %>%
  ggplot(aes(x=variable, y=.)) +
  geom_bar(stat = 'identity') +
  theme(axis.text.x = element_text(angle = 90))

Predict Social Distancing

Change point analysis

# keep only counties with full data
kreis_complete <- df_ger_socdist_scaled %>% 
  group_by(kreis) %>% 
  summarize(n = n()) %>% 
  filter(n==max(.$n)) %>% 
  .$kreis

# run changepoint analysis
df_ger_socdist_cpt_results <- df_ger_socdist_scaled %>% 
  select(kreis, socdist_single_tile) %>%
  filter(kreis %in% kreis_complete) %>% 
  split(.$kreis) %>%
  map(~ cpt.meanvar(as.vector(.$socdist_single_tile),
                    #penalty = 'Asymptotic',
                    class=TRUE,
                    param.estimates=TRUE,
                    Q=1,
                    test.stat = 'Normal'))

# calculate change point
df_ger_socdist_cpt_day <- df_ger_socdist_cpt_results %>% 
  map(cpts) %>% 
  unlist() %>% 
  as.data.frame() %>% 
  rename(cpt_day_socdist = '.') %>%
  rownames_to_column('kreis')

# calculate mean differences
df_ger_socdist_cpt_mean_diff <- df_ger_socdist_cpt_results %>% 
  map(param.est) %>% 
  map(~ .$mean) %>% 
  map(~ .[2]-.[1]) %>% 
  unlist() %>% 
  as.data.frame() %>% 
  rename(mean_diff_socdist = '.') %>%
  rownames_to_column('kreis')

# calculate varaince differences
df_ger_socdist_cpt_var_diff <- df_ger_socdist_cpt_results %>% 
  map(param.est) %>% 
  map(~ .$variance) %>% 
  map(~ .[2]-.[1]) %>% 
  unlist() %>% 
  as.data.frame() %>% 
  rename(var_diff_socdist = '.') %>%
  rownames_to_column('kreis')

# merge with county data
df_ger_cpt_socdist <- df_ger_socdist_scaled %>% 
  select(-time, -socdist_single_tile) %>%
  distinct() %>% 
  left_join(df_ger_socdist_cpt_day, by='kreis') %>%
  left_join(df_ger_socdist_cpt_mean_diff, by='kreis') %>%
  left_join(df_ger_socdist_cpt_var_diff, by='kreis') %>%
  left_join(select(df_ger_onset_prev, kreis, onset_prev), by='kreis') %>%
  left_join(select(df_ger_slope_prev, kreis, slope_prev), by='kreis') 

# handle censored data
df_ger_cpt_socdist <- df_ger_cpt_socdist %>% 
  mutate(cpt_day_socdist = ifelse(is.na(cpt_day_socdist), 
                                  as.numeric(diff(range(df_ger_socdist$date))), 
                                  cpt_day_socdist)) %>% 
  mutate(event = ifelse(cpt_day_socdist >= 60, 0, 1))

df_ger_cpt_socdist$cpt_day_socdist %>% hist()

df_ger_cpt_socdist$mean_diff_socdist %>% hist()

df_ger_cpt_socdist$var_diff_socdist %>% hist()

for(i in head(df_ger_socdist_cpt_results, 5)){
  plot(i)
}

NA

Predicting change points with time-to-event regression

# predict hazard from personality
cox_cpt_socdist <- coxph(Surv(cpt_day_socdist, event) ~ 
                           pers_o + pers_c + pers_e + pers_a + pers_n, 
                  data = df_ger_cpt_socdist)
cox_cpt_socdist %>% summary()

Call:
coxph(formula = Surv(cpt_day_socdist, event) ~ pers_o + pers_c + 
    pers_e + pers_a + pers_n, data = df_ger_cpt_socdist)

  n= 400, number of events= 400 

           coef exp(coef) se(coef)      z Pr(>|z|)  
pers_o  0.13238   1.14154  0.05555  2.383   0.0172 *
pers_c -0.04634   0.95472  0.05364 -0.864   0.3876  
pers_e -0.08928   0.91459  0.05443 -1.640   0.1009  
pers_a -0.05236   0.94898  0.04881 -1.073   0.2833  
pers_n -0.08884   0.91499  0.05376 -1.653   0.0984 .
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

       exp(coef) exp(-coef) lower .95 upper .95
pers_o    1.1415      0.876    1.0238     1.273
pers_c    0.9547      1.047    0.8594     1.061
pers_e    0.9146      1.093    0.8221     1.018
pers_a    0.9490      1.054    0.8624     1.044
pers_n    0.9150      1.093    0.8235     1.017

Concordance= 0.595  (se = 0.024 )
Likelihood ratio test= 10.81  on 5 df,   p=0.06
Wald test            = 10.76  on 5 df,   p=0.06
Score (logrank) test = 10.71  on 5 df,   p=0.06

# predict hazard from personality with controls
cox_cpt_socdist_ctrl <- coxph(Surv(cpt_day_socdist, event) ~ 
                                pers_o + pers_c + pers_e + pers_a + pers_n +
                               women + academics + cdu + afd + hospital_beds +
                               tourism_beds + gdp + manufact + airport + age +
                               popdens,
                           data = df_ger_cpt_socdist,)
cox_cpt_socdist_ctrl %>% summary()

Call:
coxph(formula = Surv(cpt_day_socdist, event) ~ pers_o + pers_c + 
    pers_e + pers_a + pers_n + women + academics + cdu + afd + 
    hospital_beds + tourism_beds + gdp + manufact + airport + 
    age + popdens, data = df_ger_cpt_socdist)

  n= 392, number of events= 392 
   (8 observations deleted due to missingness)

                   coef exp(coef)  se(coef)      z Pr(>|z|)   
pers_o         0.013954  1.014052  0.064817  0.215  0.82954   
pers_c        -0.035378  0.965241  0.061566 -0.575  0.56554   
pers_e        -0.157296  0.854451  0.061888 -2.542  0.01103 * 
pers_a        -0.008478  0.991558  0.057259 -0.148  0.88229   
pers_n        -0.043137  0.957780  0.060730 -0.710  0.47751   
women          0.084290  1.087945  0.076682  1.099  0.27167   
academics     -0.055918  0.945616  0.109764 -0.509  0.61044   
cdu           -0.135885  0.872943  0.083444 -1.628  0.10343   
afd           -0.105132  0.900205  0.090520 -1.161  0.24547   
hospital_beds -0.196237  0.821818  0.065813 -2.982  0.00287 **
tourism_beds   0.102105  1.107500  0.057115  1.788  0.07382 . 
gdp            0.135235  1.144806  0.126125  1.072  0.28362   
manufact      -0.125882  0.881719  0.093401 -1.348  0.17774   
airport       -0.074355  0.928342  0.065237 -1.140  0.25438   
age           -0.181020  0.834419  0.098405 -1.840  0.06584 . 
popdens        0.135757  1.145404  0.087636  1.549  0.12136   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

              exp(coef) exp(-coef) lower .95 upper .95
pers_o           1.0141     0.9861    0.8931    1.1514
pers_c           0.9652     1.0360    0.8555    1.0890
pers_e           0.8545     1.1703    0.7568    0.9646
pers_a           0.9916     1.0085    0.8863    1.1093
pers_n           0.9578     1.0441    0.8503    1.0788
women            1.0879     0.9192    0.9361    1.2644
academics        0.9456     1.0575    0.7626    1.1726
cdu              0.8729     1.1455    0.7412    1.0280
afd              0.9002     1.1109    0.7539    1.0750
hospital_beds    0.8218     1.2168    0.7224    0.9350
tourism_beds     1.1075     0.9029    0.9902    1.2387
gdp              1.1448     0.8735    0.8941    1.4659
manufact         0.8817     1.1341    0.7342    1.0588
airport          0.9283     1.0772    0.8169    1.0550
age              0.8344     1.1984    0.6881    1.0119
popdens          1.1454     0.8731    0.9646    1.3600

Concordance= 0.703  (se = 0.021 )
Likelihood ratio test= 64.83  on 16 df,   p=8e-08
Wald test            = 68.61  on 16 df,   p=2e-08
Score (logrank) test = 70.59  on 16 df,   p=8e-09

# predict hazard from personality with controls
cox_cpt_socdist_ctrl2 <- coxph(Surv(cpt_day_socdist, event) ~ 
                                 pers_o + pers_c + pers_e + pers_a + pers_n +
                               women + academics + cdu + afd + hospital_beds +
                               tourism_beds + gdp + manufact + airport + age +
                               popdens + onset_prev + slope_prev,
                  data = df_ger_cpt_socdist)
cox_cpt_socdist_ctrl2 %>% summary()

Call:
coxph(formula = Surv(cpt_day_socdist, event) ~ pers_o + pers_c + 
    pers_e + pers_a + pers_n + women + academics + cdu + afd + 
    hospital_beds + tourism_beds + gdp + manufact + airport + 
    age + popdens + onset_prev + slope_prev, data = df_ger_cpt_socdist)

  n= 392, number of events= 392 
   (8 observations deleted due to missingness)

                   coef exp(coef)  se(coef)      z Pr(>|z|)   
pers_o         0.011438  1.011504  0.065327  0.175  0.86100   
pers_c        -0.035315  0.965301  0.061681 -0.573  0.56695   
pers_e        -0.164389  0.848412  0.062459 -2.632  0.00849 **
pers_a        -0.001916  0.998086  0.058253 -0.033  0.97377   
pers_n        -0.026007  0.974329  0.061798 -0.421  0.67388   
women          0.095528  1.100239  0.078771  1.213  0.22524   
academics     -0.083599  0.919800  0.111618 -0.749  0.45388   
cdu           -0.148737  0.861796  0.087875 -1.693  0.09053 . 
afd           -0.092155  0.911964  0.091530 -1.007  0.31402   
hospital_beds -0.191478  0.825738  0.066089 -2.897  0.00376 **
tourism_beds   0.103020  1.108513  0.057201  1.801  0.07170 . 
gdp            0.171296  1.186842  0.127477  1.344  0.17903   
manufact      -0.151711  0.859237  0.094398 -1.607  0.10802   
airport       -0.064266  0.937756  0.065954 -0.974  0.32985   
age           -0.173100  0.841053  0.099696 -1.736  0.08251 . 
popdens        0.121096  1.128733  0.087691  1.381  0.16730   
onset_prev    -0.008623  0.991414  0.005259 -1.640  0.10110   
slope_prev     0.259984  1.296909  2.375477  0.109  0.91285   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

              exp(coef) exp(-coef) lower .95 upper .95
pers_o           1.0115     0.9886   0.88994    1.1497
pers_c           0.9653     1.0359   0.85538    1.0893
pers_e           0.8484     1.1787   0.75066    0.9589
pers_a           0.9981     1.0019   0.89040    1.1188
pers_n           0.9743     1.0263   0.86318    1.0998
women            1.1002     0.9089   0.94284    1.2839
academics        0.9198     1.0872   0.73907    1.1447
cdu              0.8618     1.1604   0.72545    1.0238
afd              0.9120     1.0965   0.76220    1.0912
hospital_beds    0.8257     1.2110   0.72542    0.9399
tourism_beds     1.1085     0.9021   0.99095    1.2400
gdp              1.1868     0.8426   0.92445    1.5237
manufact         0.8592     1.1638   0.71410    1.0339
airport          0.9378     1.0664   0.82404    1.0672
age              0.8411     1.1890   0.69177    1.0226
popdens          1.1287     0.8859   0.95049    1.3404
onset_prev       0.9914     1.0087   0.98125    1.0017
slope_prev       1.2969     0.7711   0.01233  136.4328

Concordance= 0.704  (se = 0.021 )
Likelihood ratio test= 68.07  on 18 df,   p=1e-07
Wald test            = 71.78  on 18 df,   p=2e-08
Score (logrank) test = 73.73  on 18 df,   p=1e-08

Linear models predicting mean differences

lm_meandiff_socdist <- lm(mean_diff_socdist ~ 
                            pers_o + pers_c + pers_e + pers_a + pers_n, 
                         data = df_ger_cpt_socdist)
lm_meandiff_socdist %>% summary()

Call:
lm(formula = mean_diff_socdist ~ pers_o + pers_c + pers_e + pers_a + 
    pers_n, data = df_ger_cpt_socdist)

Residuals:
      Min        1Q    Median        3Q       Max 
-0.042484 -0.009668 -0.000179  0.008346  0.062165 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  0.0949484  0.0007323 129.654  < 2e-16 ***
pers_o       0.0056844  0.0007839   7.251 2.21e-12 ***
pers_c      -0.0008564  0.0008316  -1.030   0.3037    
pers_e       0.0008648  0.0008356   1.035   0.3014    
pers_a      -0.0002079  0.0008311  -0.250   0.8026    
pers_n      -0.0020169  0.0008627  -2.338   0.0199 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.01465 on 394 degrees of freedom
Multiple R-squared:  0.1631,    Adjusted R-squared:  0.1525 
F-statistic: 15.36 on 5 and 394 DF,  p-value: 8.43e-14

lm_meandiff_socdist_ctrl <- lm(mean_diff_socdist ~ 
                                 pers_o + pers_c + pers_e + pers_a + pers_n +
                               women + academics + cdu + afd + hospital_beds +
                               tourism_beds + gdp + manufact + airport + age +
                               popdens,
                           data = df_ger_cpt_socdist,)
lm_meandiff_socdist_ctrl %>% summary()

Call:
lm(formula = mean_diff_socdist ~ pers_o + pers_c + pers_e + pers_a + 
    pers_n + women + academics + cdu + afd + hospital_beds + 
    tourism_beds + gdp + manufact + airport + age + popdens, 
    data = df_ger_cpt_socdist)

Residuals:
      Min        1Q    Median        3Q       Max 
-0.045968 -0.007500 -0.000700  0.006643  0.039513 

Coefficients:
                Estimate Std. Error t value Pr(>|t|)    
(Intercept)    9.485e-02  5.814e-04 163.137  < 2e-16 ***
pers_o         1.708e-03  7.178e-04   2.379 0.017838 *  
pers_c         7.003e-04  6.958e-04   1.007 0.314809    
pers_e        -8.269e-04  6.850e-04  -1.207 0.228176    
pers_a        -3.770e-04  7.047e-04  -0.535 0.592923    
pers_n        -2.724e-05  7.107e-04  -0.038 0.969449    
women          2.286e-03  8.574e-04   2.666 0.008013 ** 
academics      1.987e-03  1.157e-03   1.716 0.086908 .  
cdu            2.873e-03  9.205e-04   3.121 0.001942 ** 
afd           -1.348e-03  9.043e-04  -1.490 0.137006    
hospital_beds -2.819e-03  7.821e-04  -3.604 0.000355 ***
tourism_beds   4.121e-05  6.670e-04   0.062 0.950768    
gdp            5.408e-03  1.343e-03   4.028 6.82e-05 ***
manufact      -1.267e-03  1.079e-03  -1.175 0.240741    
airport       -1.163e-03  7.317e-04  -1.589 0.112832    
age           -4.297e-03  1.118e-03  -3.842 0.000143 ***
popdens        7.862e-04  9.986e-04   0.787 0.431608    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.0115 on 375 degrees of freedom
  (8 observations deleted due to missingness)
Multiple R-squared:  0.5001,    Adjusted R-squared:  0.4787 
F-statistic: 23.44 on 16 and 375 DF,  p-value: < 2.2e-16

lm_meandiff_socdist_ctrl2 <- lm(mean_diff_socdist ~ 
                                  pers_o + pers_c + pers_e + pers_a + pers_n +
                               women + academics + cdu + afd + hospital_beds +
                               tourism_beds + gdp + manufact + airport + age +
                               popdens + onset_prev + slope_prev,
                  data = df_ger_cpt_socdist)
lm_meandiff_socdist_ctrl2 %>% summary()

Call:
lm(formula = mean_diff_socdist ~ pers_o + pers_c + pers_e + pers_a + 
    pers_n + women + academics + cdu + afd + hospital_beds + 
    tourism_beds + gdp + manufact + airport + age + popdens + 
    onset_prev + slope_prev, data = df_ger_cpt_socdist)

Residuals:
      Min        1Q    Median        3Q       Max 
-0.044599 -0.007374 -0.000808  0.006306  0.040989 

Coefficients:
                Estimate Std. Error t value Pr(>|t|)    
(Intercept)    9.476e-02  3.013e-03  31.446  < 2e-16 ***
pers_o         1.742e-03  7.021e-04   2.481 0.013540 *  
pers_c         7.559e-04  6.818e-04   1.109 0.268258    
pers_e        -9.932e-04  6.721e-04  -1.478 0.140321    
pers_a        -4.823e-04  6.900e-04  -0.699 0.485022    
pers_n        -2.376e-05  7.018e-04  -0.034 0.973012    
women          1.794e-03  8.465e-04   2.119 0.034720 *  
academics      1.794e-03  1.142e-03   1.571 0.117084    
cdu            1.747e-03  9.532e-04   1.833 0.067665 .  
afd           -1.604e-03  8.864e-04  -1.809 0.071255 .  
hospital_beds -3.013e-03  7.706e-04  -3.910 0.000109 ***
tourism_beds   3.179e-04  6.567e-04   0.484 0.628597    
gdp            5.439e-03  1.317e-03   4.131 4.47e-05 ***
manufact      -1.618e-03  1.060e-03  -1.526 0.127771    
airport       -1.478e-03  7.197e-04  -2.054 0.040649 *  
age           -4.000e-03  1.099e-03  -3.640 0.000312 ***
popdens        7.496e-04  9.795e-04   0.765 0.444575    
onset_prev    -5.130e-05  5.793e-05  -0.885 0.376476    
slope_prev     1.077e-01  2.593e-02   4.152 4.08e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.01125 on 373 degrees of freedom
  (8 observations deleted due to missingness)
Multiple R-squared:  0.5243,    Adjusted R-squared:  0.5014 
F-statistic: 22.84 on 18 and 373 DF,  p-value: < 2.2e-16

CRF predicting mean difference

ctrls <- cforest_unbiased(ntree=500, mtry=5)

crf_meandiff_socdist <- cforest(mean_diff_socdist ~ 
                                  pers_o + pers_c + pers_e + pers_a + pers_n +
                               women + academics + cdu + afd + hospital_beds +
                               tourism_beds + gdp + manufact + airport + age +
                               popdens + onset_prev + slope_prev,
                              data = df_ger_cpt_socdist %>% drop_na(),
                         controls = ctrls)

crf_meandiff_socdist_varimp <- varimp(crf_meandiff_socdist, nperm = 1)
crf_meandiff_socdist_varimp_cond <- varimp(crf_meandiff_socdist, conditional = T, nperm = 1)

crf_meandiff_socdist_varimp %>% as.data.frame() %>% 
  rownames_to_column('variable') %>%
  ggplot(aes(x=variable, y=.)) +
  geom_bar(stat = 'identity') +
  theme(axis.text.x = element_text(angle = 90))

crf_meandiff_socdist_varimp_cond %>% as.data.frame() %>% 
  rownames_to_column('variable') %>%
  ggplot(aes(x=variable, y=.)) +
  geom_bar(stat = 'identity') +
  theme(axis.text.x = element_text(angle = 90))

Export data

us_list_results <- list(cox_onset_prev, cox_onset_prev_ctrl, 
     lm_slope_prev, lm_slope_prev_coef, 
     cox_cpt_socdist, cox_cpt_socdist_ctrl, cox_cpt_socdist_ctrl2,
     lm_meandiff_socdist, lm_meandiff_socdist_ctrl, lm_meandiff_socdist_ctrl2)

results_names <- list('cox_onset_prev', 'cox_onset_prev_ctrl', 
     'lm_slope_prev', 'lm_slope_prev_coef', 
     'cox_cpt_socdist', 'cox_cpt_socdist_ctrl', 'cox_cpt_socdist_ctrl2',
     'lm_meandiff_socdist', 'lm_meandiff_socdist_ctrl', 'lm_meandiff_socdist_ctrl2')

names(us_list_results) <- results_names

save(us_list_results, file="us_list_results.RData")

write_csv(df_us_slope_prev, 'df_us_slope_prev.csv')
write_csv(df_us_cpt_socdist, 'df_us_cpt_socdist.csv')